
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    
    <title>odespy &mdash; Odespy API 0.2 documentation</title>
    
    <link rel="stylesheet" href="_static/pyramid.css" type="text/css" />
    <link rel="stylesheet" href="_static/pygments.css" type="text/css" />
    
    <script type="text/javascript">
      var DOCUMENTATION_OPTIONS = {
        URL_ROOT:    '',
        VERSION:     '0.2',
        COLLAPSE_INDEX: false,
        FILE_SUFFIX: '.html',
        HAS_SOURCE:  true
      };
    </script>
    <script type="text/javascript" src="_static/jquery.js"></script>
    <script type="text/javascript" src="_static/underscore.js"></script>
    <script type="text/javascript" src="_static/doctools.js"></script>
    <script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
    <link rel="top" title="Odespy API 0.2 documentation" href="index.html" />
    <link rel="next" title="odespy.rkc" href="rkc.html" />
    <link rel="prev" title="Odespy API" href="index.html" />
<link rel="stylesheet" href="http://fonts.googleapis.com/css?family=Neuton&amp;subset=latin" type="text/css" media="screen" charset="utf-8" />
<link rel="stylesheet" href="http://fonts.googleapis.com/css?family=Nobile:regular,italic,bold,bolditalic&amp;subset=latin" type="text/css" media="screen" charset="utf-8" />
<!--[if lte IE 6]>
<link rel="stylesheet" href="_static/ie6.css" type="text/css" media="screen" charset="utf-8" />
<![endif]-->

  </head>
  <body>

    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             accesskey="I">index</a></li>
        <li class="right" >
          <a href="py-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li class="right" >
          <a href="np-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li class="right" >
          <a href="rkc.html" title="odespy.rkc"
             accesskey="N">next</a> |</li>
        <li class="right" >
          <a href="index.html" title="Odespy API"
             accesskey="P">previous</a> |</li>
        <li><a href="index.html">Odespy API 0.2 documentation</a> &raquo;</li> 
      </ul>
    </div>  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body">
            
  <div class="section" id="module-odespy">
<span id="odespy"></span><h1><a class="reference internal" href="#module-odespy" title="odespy"><tt class="xref py py-mod docutils literal"><span class="pre">odespy</span></tt></a><a class="headerlink" href="#module-odespy" title="Permalink to this headline">¶</a></h1>
<p>The <tt class="docutils literal"><span class="pre">odespy</span></tt> package contains tools for solving ordinary
differential equations (ODEs). The user specifies the problem through
high-level Python code. Both scalar ODEs and systems of ODEs are
supported.  A wide range of numerical methods for ODEs are offered:</p>
<table border="1" class="docutils">
<colgroup>
<col width="24%" />
<col width="76%" />
</colgroup>
<thead valign="bottom">
<tr class="row-odd"><th class="head">Classname</th>
<th class="head">Short description</th>
</tr>
</thead>
<tbody valign="top">
<tr class="row-even"><td>AdamsBashMoulton2</td>
<td>Explicit 2nd-order Adams-Bashforth-Moulton method</td>
</tr>
<tr class="row-odd"><td>AdamsBashMoulton3</td>
<td>Explicit 3rd-order Adams-Bashforth-Moulton method</td>
</tr>
<tr class="row-even"><td>AdamsBashforth2</td>
<td>Explicit 2nd-order Adams-Bashforth method</td>
</tr>
<tr class="row-odd"><td>AdamsBashforth3</td>
<td>Explicit 3rd-order Adams-Bashforth method</td>
</tr>
<tr class="row-even"><td>AdamsBashforth4</td>
<td>Explicit 4th-order Adams-Bashforth method</td>
</tr>
<tr class="row-odd"><td>AdaptiveResidual</td>
<td>Very simple adaptive strategy based on the residual</td>
</tr>
<tr class="row-even"><td>Backward2Step</td>
<td>Implicit 2nd-order Backward Euler method</td>
</tr>
<tr class="row-odd"><td>BackwardEuler</td>
<td>Implicit 1st-order Backward Euler method</td>
</tr>
<tr class="row-even"><td>BogackiShampine</td>
<td>Adaptive Bogacki-Shampine RK method of order (3, 2)</td>
</tr>
<tr class="row-odd"><td>CashKarp</td>
<td>Adaptive Cash-Karp RK method of order (5, 4)</td>
</tr>
<tr class="row-even"><td>Dop853</td>
<td>Adaptive Dormand &amp; Prince method of order 8(5,3) (scipy)</td>
</tr>
<tr class="row-odd"><td>Dopri5</td>
<td>Dormand &amp; Prince method of order 5(4) (scipy)</td>
</tr>
<tr class="row-even"><td>DormandPrince</td>
<td>Dormand &amp; Prince RK method of order (5, 4)</td>
</tr>
<tr class="row-odd"><td>Euler</td>
<td>The simple explicit (forward) Euler scheme</td>
</tr>
<tr class="row-even"><td>Fehlberg</td>
<td>Adaptive Runge-Kutta-Fehlberg (4,5) method</td>
</tr>
<tr class="row-odd"><td>ForwardEuler</td>
<td>The simple explicit (forward) Euler scheme</td>
</tr>
<tr class="row-even"><td>Heun</td>
<td>Heun&#8217;s explicit method (similar to RK2)</td>
</tr>
<tr class="row-odd"><td>Leapfrog</td>
<td>Standard explicit Leapfrog scheme</td>
</tr>
<tr class="row-even"><td>LeapfrogFiltered</td>
<td>Filtered Leapfrog scheme</td>
</tr>
<tr class="row-odd"><td>Lsoda</td>
<td>LSODA solver with stiff-nonstiff auto shift</td>
</tr>
<tr class="row-even"><td>Lsodar</td>
<td>LSODAR method with stiff-nonstiff auto shift</td>
</tr>
<tr class="row-odd"><td>Lsode</td>
<td>LSODE solver for a stiff or nonstiff system</td>
</tr>
<tr class="row-even"><td>Lsodes</td>
<td>LSODES solver for sparse Jacobians</td>
</tr>
<tr class="row-odd"><td>Lsodi</td>
<td>LSODI solver for linearly implicit systems</td>
</tr>
<tr class="row-even"><td>Lsodis</td>
<td>LSODIS solver for linearly implicit sparse systems</td>
</tr>
<tr class="row-odd"><td>Lsoibt</td>
<td>LSOIBIT solver for linearly implicit block tridiag systems</td>
</tr>
<tr class="row-even"><td>MidpointImplicit</td>
<td>Implicit 2nd-order Midpoint method</td>
</tr>
<tr class="row-odd"><td>MidpointIter</td>
<td>Explicit 2nd-order iterated Midpoint method</td>
</tr>
<tr class="row-even"><td>RK2</td>
<td>Explicit 2nd-order Runge-Kutta method</td>
</tr>
<tr class="row-odd"><td>RK3</td>
<td>Explicit 3rd-order Runge-Kutta method</td>
</tr>
<tr class="row-even"><td>RK34</td>
<td>Adaptive 4th-order Runge-Kutta method</td>
</tr>
<tr class="row-odd"><td>RK4</td>
<td>Explicit 4th-order Runge-Kutta method</td>
</tr>
<tr class="row-even"><td>RKC</td>
<td>Explicit 2nd-order Runge-Kutta-Chebyshev method (rkc.f)</td>
</tr>
<tr class="row-odd"><td>RKF45</td>
<td>Adaptive Runge-Kutta-Fehlberg (4,5) method (rkf45.f)</td>
</tr>
<tr class="row-even"><td>RKFehlberg</td>
<td>Adaptive Runge-Kutta-Fehlberg (4,5) method</td>
</tr>
<tr class="row-odd"><td>RungeKutta1</td>
<td>Explicit 1st-order Runge-Kutta method</td>
</tr>
<tr class="row-even"><td>RungeKutta2</td>
<td>Explicit 2nd-order Runge-Kutta method</td>
</tr>
<tr class="row-odd"><td>RungeKutta3</td>
<td>Explicit 3rd-order Runge-Kutta method</td>
</tr>
<tr class="row-even"><td>RungeKutta4</td>
<td>Explicit 4th-order Runge-Kutta method</td>
</tr>
<tr class="row-odd"><td>ThetaRule</td>
<td>Unified Forward/Backward Euler and Midpoint methods</td>
</tr>
<tr class="row-even"><td>Trapezoidal</td>
<td>Heun&#8217;s explicit method (similar to RK2)</td>
</tr>
<tr class="row-odd"><td>Vode</td>
<td>Adams/BDF Vode adaptive method (vode.f wrapper)</td>
</tr>
<tr class="row-even"><td>lsoda_scipy</td>
<td>Wrapper of lsoda (scipy.integrate.odeint)</td>
</tr>
<tr class="row-odd"><td>odefun_sympy</td>
<td>Very accurate high order Taylor method (from SymPy)</td>
</tr>
<tr class="row-even"><td>odelab</td>
<td>interface to all solvers in odelab</td>
</tr>
</tbody>
</table>
<div class="section" id="basic-usage">
<h2>Basic Usage<a class="headerlink" href="#basic-usage" title="Permalink to this headline">¶</a></h2>
<p>This section explains how to use Odespy. The general principles and
program steps are first explained and then followed by a series of
examples with progressive complexity with respect to Python constructs
and numerical methods.</p>
<div class="section" id="overview">
<h3>Overview<a class="headerlink" href="#overview" title="Permalink to this headline">¶</a></h3>
<p>A code using Odespy to solve ODEs consists of six steps. These are
outlined in generic form below.</p>
<div class="section" id="step-1">
<h4>Step 1<a class="headerlink" href="#step-1" title="Permalink to this headline">¶</a></h4>
<p><em>Write the ODE problem in generic form</em> <span class="math">\(u' = f(u, t)\)</span>,
where <span class="math">\(u(t)\)</span> is the unknown function to be solved for, or a vector
of unknown functions of time in case of a system of ODEs.</p>
</div>
<div class="section" id="step-2">
<h4>Step 2<a class="headerlink" href="#step-2" title="Permalink to this headline">¶</a></h4>
<p><em>Implement the right-hand side function</em> <span class="math">\(f(u, t)\)</span> as a Python function
<tt class="docutils literal"><span class="pre">f(u,</span> <span class="pre">t)</span></tt>.  The argument <tt class="docutils literal"><span class="pre">u</span></tt> is either a <tt class="docutils literal"><span class="pre">float</span></tt> object, in case of a
scalar ODE, or a <tt class="docutils literal"><span class="pre">numpy</span></tt> array object, in case of a system of ODEs.
Some solvers in this package also allow implementation of <span class="math">\(f\)</span> in
FORTRAN for increased efficiency.</p>
</div>
<div class="section" id="step-3">
<h4>Step 3<a class="headerlink" href="#step-3" title="Permalink to this headline">¶</a></h4>
<p><em>Create a solver object</em></p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">solver</span> <span class="o">=</span> <span class="n">classname</span><span class="p">(</span><span class="n">f</span><span class="p">)</span>
</pre></div>
</div>
<p>where <tt class="docutils literal"><span class="pre">classname</span></tt> is the name of a class in this package implementing
the desired numerical method.</p>
<p>Many solver classes has a range of parameters that the user can set to
control various parts of the solution process. The parameters are
documented in the doc string of the class (<tt class="docutils literal"><span class="pre">pydoc</span> <span class="pre">classname</span></tt> will list
the documentation in a terminal window). One can either specify parameters
at construction time, via extra keyword arguments to the constructor,</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">solver</span> <span class="o">=</span> <span class="n">classname</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="n">prm1</span><span class="o">=</span><span class="n">value1</span><span class="p">,</span> <span class="n">prm2</span><span class="o">=</span><span class="n">value2</span><span class="p">,</span> <span class="o">...</span><span class="p">)</span>
</pre></div>
</div>
<p>or at any time using the <tt class="docutils literal"><span class="pre">set</span></tt> method:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">solver</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="n">prm1</span><span class="o">=</span><span class="n">value1</span><span class="p">,</span> <span class="n">prm2</span><span class="o">=</span><span class="n">value2</span><span class="p">,</span> <span class="n">prm3</span><span class="o">=</span><span class="n">value3</span><span class="p">)</span>
<span class="o">...</span>
<span class="n">solver</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="n">prm4</span><span class="o">=</span><span class="n">value4</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="section" id="step-4">
<h4>Step 4<a class="headerlink" href="#step-4" title="Permalink to this headline">¶</a></h4>
<p><em>Set the initial condition</em> <span class="math">\(u(0)=U_0\)</span>,</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">solver</span><span class="o">.</span><span class="n">set_initial_condition</span><span class="p">(</span><span class="n">U0</span><span class="p">)</span>
</pre></div>
</div>
<p>where <tt class="docutils literal"><span class="pre">U0</span></tt> is either a number, for a scalar ODE, or a sequence (list, tuple,
<tt class="docutils literal"><span class="pre">numpy</span></tt> array), for a system of ODEs.</p>
</div>
<div class="section" id="step-5">
<h4>Step 5<a class="headerlink" href="#step-5" title="Permalink to this headline">¶</a></h4>
<p><em>Solve the ODE problem</em>, which means to compute <span class="math">\(u(t)\)</span> at
some discrete user-specified time points <span class="math">\(t_1, t_2, \ldots, t_N\)</span>.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">T</span> <span class="o">=</span> <span class="o">...</span>  <span class="c"># end time</span>
<span class="n">time_points</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">T</span><span class="p">,</span> <span class="n">N</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>
<span class="n">u</span><span class="p">,</span> <span class="n">t</span> <span class="o">=</span> <span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">time_points</span><span class="p">)</span>
</pre></div>
</div>
<p>In case of a scalar ODE, the returned solution <tt class="docutils literal"><span class="pre">u</span></tt> is a one-dimensional
<tt class="docutils literal"><span class="pre">numpy</span></tt> array where <tt class="docutils literal"><span class="pre">u[i]</span></tt> holds the solution at time point <tt class="docutils literal"><span class="pre">t[i]</span></tt>.
For a system of ODEs, the returned <tt class="docutils literal"><span class="pre">u</span></tt> is a two-dimensional <tt class="docutils literal"><span class="pre">numpy</span></tt>
array where <tt class="docutils literal"><span class="pre">u[i,j]</span></tt> holds the solution of the $j$-th unknown
function at the $i$-th time point <tt class="docutils literal"><span class="pre">t[i]</span></tt> (<span class="math">\(u_j(t_i)\)</span> in mathematics
notation).</p>
<p>By giving the parameter <tt class="docutils literal"><span class="pre">disk_storage=True</span></tt> to the solver&#8217;s constructor,
the returned <tt class="docutils literal"><span class="pre">u</span></tt> array is memory mapped (i.e., of type <tt class="docutils literal"><span class="pre">numpy.memmap</span></tt>)
such that all the data are stored on file, but parts of the array can
be efficiently accessed.</p>
<p>The <tt class="docutils literal"><span class="pre">time_points</span></tt> array specifies the time points where we want the
solution to be computed. The returned array <tt class="docutils literal"><span class="pre">t</span></tt> is the same as
<tt class="docutils literal"><span class="pre">time_points</span></tt>.  The simplest numerical methods in the Odespy
package apply the <tt class="docutils literal"><span class="pre">time_points</span></tt> array directly for the time stepping.
That is, the time steps used are given by</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">time_points</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">-</span> <span class="n">time_points</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>  <span class="c">#  i=0,1,...,len(time_points)-1</span>
</pre></div>
</div>
<p>The adaptive schemes typically compute between each time point in
the <tt class="docutils literal"><span class="pre">time_points</span></tt> array, making this array a specification where
values of the unknowns are desired.</p>
<p>The <tt class="docutils literal"><span class="pre">solve</span></tt> method in solver classes also allows a second argument,
<tt class="docutils literal"><span class="pre">terminate</span></tt>, which is a user-implemented Python function specifying
when the solution process is to be terminated. For example,
terminating when the solution reaches an asymptotic (known) value
<tt class="docutils literal"><span class="pre">a</span></tt> can be done by</p>
<div class="highlight-python"><div class="highlight"><pre><span class="k">def</span> <span class="nf">terminate</span><span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">step_no</span><span class="p">):</span>
    <span class="c"># u and t are arrays. Most recent solution is u[step_no].</span>
    <span class="n">tolerance</span> <span class="o">=</span> <span class="mf">1E-6</span>
    <span class="k">return</span> <span class="nb">abs</span><span class="p">(</span><span class="n">u</span><span class="p">[</span><span class="n">step_no</span><span class="p">]</span> <span class="o">-</span> <span class="n">a</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">tolerance</span>

<span class="n">u</span><span class="p">,</span> <span class="n">t</span> <span class="o">=</span> <span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">time_points</span><span class="p">,</span> <span class="n">terminate</span><span class="p">)</span>
</pre></div>
</div>
<p>The arguments transferred to the <tt class="docutils literal"><span class="pre">terminate</span></tt> function are the
solution array <tt class="docutils literal"><span class="pre">u</span></tt>, the corresponding time points <tt class="docutils literal"><span class="pre">t</span></tt>, and
an integer <tt class="docutils literal"><span class="pre">step_no</span></tt> reflecting the most recently computed <tt class="docutils literal"><span class="pre">u</span></tt>
value. That is, <tt class="docutils literal"><span class="pre">u[step_no]</span></tt> is most recently computed value of <span class="math">\(u\)</span>.
(The array data <tt class="docutils literal"><span class="pre">u[step_no+1:]</span></tt> will typically be zero as these
are uncomputed future values.)</p>
</div>
<div class="section" id="step-6">
<h4>Step 6<a class="headerlink" href="#step-6" title="Permalink to this headline">¶</a></h4>
<p><em>Extract solution components</em> for plotting and further analysis.
Since the <tt class="docutils literal"><span class="pre">u</span></tt> array returned from <tt class="docutils literal"><span class="pre">solver.solve</span></tt> stores all unknown
functions at all discrete time levels, one usually wants to extract
individual unknowns as one-dimensional arrays. Here is an example
where unknown number <span class="math">\(0\)</span> and <span class="math">\(k\)</span> are extracted in individual arrays
and plotted:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">u_0</span> <span class="o">=</span> <span class="n">u</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span>
<span class="n">u_k</span> <span class="o">=</span> <span class="n">u</span><span class="p">[:,</span><span class="n">k</span><span class="p">]</span>

<span class="kn">from</span> <span class="nn">matplotlib.pyplot</span> <span class="kn">import</span> <span class="n">plot</span><span class="p">,</span> <span class="n">show</span>
<span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span> <span class="n">u_0</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">u_k</span><span class="p">)</span>
<span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
</div>
</div>
<div class="section" id="first-example-logistic-growth">
<span id="ode-sec-exgr"></span><h3>First Example: Logistic Growth<a class="headerlink" href="#first-example-logistic-growth" title="Permalink to this headline">¶</a></h3>
<p>Our first example concerns the simple scalar ODE problem</p>
<div class="math">
\[\frac{du}{dt}=au\left(1-\frac{u}{R}\right),\quad u(0)=A,\]</div>
<p>where <span class="math">\(A&gt;0\)</span>, <span class="math">\(a&gt;0\)</span>, and <span class="math">\(R&gt;0\)</span> are known constants. This is
a common model for population dynamics in ecology where <span class="math">\(u\)</span> is the
number of individuals, <span class="math">\(a\)</span> the initial growth rate, <span class="math">\(R\)</span> is the
maximum number of individuals that the environment allows (the so-called
<em>carrying capacity</em> of the environment).</p>
<p>Using a standard
Runge-Kutta method of order four, the code for solving the problem in
the time interval <span class="math">\([0,10]\)</span> with <span class="math">\(N=30\)</span> time steps, looks like this
(program <a class="reference external" href="https://github.com/hplgit/odespy/blob/master/doc/src/odespy/src-odespy/logistic1.py">logistic1.py</a>):</p>
<div class="highlight-python"><div class="highlight"><pre><span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
    <span class="k">return</span> <span class="n">a</span><span class="o">*</span><span class="n">u</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u</span><span class="o">/</span><span class="n">R</span><span class="p">)</span>

<span class="n">a</span> <span class="o">=</span> <span class="mi">2</span>
<span class="n">R</span> <span class="o">=</span> <span class="mf">1E+5</span>
<span class="n">A</span> <span class="o">=</span> <span class="mi">1</span>

<span class="kn">import</span> <span class="nn">odespy</span>
<span class="n">solver</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">RK4</span><span class="p">(</span><span class="n">f</span><span class="p">)</span>
<span class="n">solver</span><span class="o">.</span><span class="n">set_initial_condition</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>

<span class="kn">from</span> <span class="nn">numpy</span> <span class="kn">import</span> <span class="n">linspace</span><span class="p">,</span> <span class="n">exp</span>
<span class="n">T</span> <span class="o">=</span> <span class="mi">10</span>  <span class="c"># end of simulation</span>
<span class="n">N</span> <span class="o">=</span> <span class="mi">30</span>  <span class="c"># no of time steps</span>
<span class="n">time_points</span> <span class="o">=</span> <span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">T</span><span class="p">,</span> <span class="n">N</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>
<span class="n">u</span><span class="p">,</span> <span class="n">t</span> <span class="o">=</span> <span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">time_points</span><span class="p">)</span>
</pre></div>
</div>
<p>With the <tt class="docutils literal"><span class="pre">RK4</span></tt> method and other non-adaptive methods
the time steps are dictated by the <tt class="docutils literal"><span class="pre">time_points</span></tt> array.
A constant time step of size is implied in the present example.
Running an alternative numerical method just means replacing <tt class="docutils literal"><span class="pre">RK4</span></tt> by, e.g.,
<tt class="docutils literal"><span class="pre">RK2</span></tt>, <tt class="docutils literal"><span class="pre">ForwardEuler</span></tt>, <tt class="docutils literal"><span class="pre">BackwardEuler</span></tt>, <tt class="docutils literal"><span class="pre">AdamsBashforth2</span></tt>, etc.</p>
<p>We can easily plot the numerical solution and compare with the exact
solution (which is known for this equation):</p>
<div class="highlight-python"><div class="highlight"><pre><span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="n">t</span><span class="p">):</span>
    <span class="k">return</span> <span class="n">R</span><span class="o">*</span><span class="n">A</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="o">*</span><span class="n">t</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">R</span> <span class="o">+</span> <span class="n">A</span><span class="o">*</span><span class="p">(</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="o">*</span><span class="n">t</span><span class="p">)</span> <span class="o">-</span> <span class="mi">1</span><span class="p">))</span>

<span class="kn">from</span> <span class="nn">matplotlib.pyplot</span> <span class="kn">import</span> <span class="o">*</span>

<span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="s">&#39;r-&#39;</span><span class="p">,</span>
     <span class="n">t</span><span class="p">,</span> <span class="n">u_exact</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="s">&#39;bo&#39;</span><span class="p">)</span>
<span class="n">savefig</span><span class="p">(</span><span class="s">&#39;tmp.png&#39;</span><span class="p">)</span>
<span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
<div class="figure">
<img alt="_images/logistic1.png" src="_images/logistic1.png" style="width: 400px;" />
<p class="caption">Solution of the logistic equation with the 4-th order Runge-Kutta method (solid line) and comparison with the exact solution (dots)</p>
</div>
<p>All the examples in this tutorial are found in the GitHub directory
<a class="reference external" href="https://github.com/hplgit/odespy/tree/master/doc/src/odespy/src-odespy/">https://github.com/hplgit/odespy/tree/master/doc/src/odespy/src-odespy/</a>.
If you download the tarball or clone the GitHub repository, the examples
reside in the directory <tt class="docutils literal"><span class="pre">doc/src/odespy/src-odespy</span></tt>.</p>
</div>
<div class="section" id="parameters-in-the-right-hand-side-function">
<span id="ode-sec-exgr-farg"></span><h3>Parameters in the Right-Hand Side Function<a class="headerlink" href="#parameters-in-the-right-hand-side-function" title="Permalink to this headline">¶</a></h3>
<p>The right-hand side function and all physical parameters are often
lumped together in a class, for instance,</p>
<div class="highlight-python"><div class="highlight"><pre><span class="k">class</span> <span class="nc">Logistic</span><span class="p">:</span>
    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">R</span><span class="p">,</span> <span class="n">A</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">a</span> <span class="o">=</span> <span class="n">a</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">R</span> <span class="o">=</span> <span class="n">R</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">A</span> <span class="o">=</span> <span class="n">A</span>

    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">a</span><span class="p">,</span> <span class="n">R</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">R</span>  <span class="c"># short form</span>
        <span class="k">return</span> <span class="n">a</span><span class="o">*</span><span class="n">u</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u</span><span class="o">/</span><span class="n">R</span><span class="p">)</span>

    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">a</span><span class="p">,</span> <span class="n">R</span><span class="p">,</span> <span class="n">A</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">R</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">A</span>  <span class="c"># short form</span>
        <span class="k">return</span> <span class="n">R</span><span class="o">*</span><span class="n">A</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="o">*</span><span class="n">t</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">R</span> <span class="o">+</span> <span class="n">A</span><span class="o">*</span><span class="p">(</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="o">*</span><span class="n">t</span><span class="p">)</span> <span class="o">-</span> <span class="mi">1</span><span class="p">))</span>
</pre></div>
</div>
<p>Note that introducing local variables like <tt class="docutils literal"><span class="pre">a</span></tt> and <tt class="docutils literal"><span class="pre">R</span></tt>, instead of
using <tt class="docutils literal"><span class="pre">self.a</span></tt> and <tt class="docutils literal"><span class="pre">self.A</span></tt>, makes the code closer to the mathematics.
This can be convenient when proof reading the implementation of
complicated ODEs.</p>
<p>The numerical solution is computed by</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">odespy</span>
<span class="n">problem</span> <span class="o">=</span> <span class="n">Logistic</span><span class="p">(</span><span class="n">a</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span> <span class="n">R</span><span class="o">=</span><span class="mf">1E+5</span><span class="p">,</span> <span class="n">A</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
<span class="n">solver</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">RK4</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">f</span><span class="p">)</span>
<span class="n">solver</span><span class="o">.</span><span class="n">set_initial_condition</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">A</span><span class="p">)</span>

<span class="n">T</span> <span class="o">=</span> <span class="mi">10</span>  <span class="c"># end of simulation</span>
<span class="n">N</span> <span class="o">=</span> <span class="mi">30</span>  <span class="c"># no of time steps</span>
<span class="n">time_points</span> <span class="o">=</span> <span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">T</span><span class="p">,</span> <span class="n">N</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>
<span class="n">u</span><span class="p">,</span> <span class="n">t</span> <span class="o">=</span> <span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">time_points</span><span class="p">)</span>
</pre></div>
</div>
<p>The complete program is available in the file program <a class="reference external" href="https://github.com/hplgit/odespy/blob/master/doc/src/odespy/src-odespy/logistic2.py">logistic2.py</a>.</p>
<p>Instead of having the problem parameters <tt class="docutils literal"><span class="pre">a</span></tt> and <tt class="docutils literal"><span class="pre">R</span></tt> in the ODE as
global variables or in a class, we may include them as extra arguments
to <tt class="docutils literal"><span class="pre">f</span></tt>, either as positional arguments or as keyword
arguments. Positional arguments can be sent to <tt class="docutils literal"><span class="pre">f</span></tt> via the constructor
argument <tt class="docutils literal"><span class="pre">f_args</span></tt> (a list/tuple of variables), while a dictionary
<tt class="docutils literal"><span class="pre">f_kwargs</span></tt> is used to transfer keyword arguments to <tt class="docutils literal"><span class="pre">f</span></tt> via the
constructor. Here is an example on using keyword arguments:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">R</span><span class="o">=</span><span class="mi">1</span><span class="p">):</span>
    <span class="k">return</span> <span class="n">a</span><span class="o">*</span><span class="n">u</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u</span><span class="o">/</span><span class="n">R</span><span class="p">)</span>

<span class="n">A</span> <span class="o">=</span> <span class="mi">1</span>

<span class="kn">import</span> <span class="nn">odespy</span>
<span class="n">solver</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">RK4</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="n">f_kwargs</span><span class="o">=</span><span class="nb">dict</span><span class="p">(</span><span class="n">a</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span> <span class="n">R</span><span class="o">=</span><span class="mf">1E+5</span><span class="p">))</span>
</pre></div>
</div>
<p>In general, a mix
of positional and keyword arguments can be used in <tt class="docutils literal"><span class="pre">f</span></tt>:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">arg1</span><span class="p">,</span> <span class="n">arg2</span><span class="p">,</span> <span class="n">arg3</span><span class="p">,</span> <span class="o">...</span><span class="p">,</span> <span class="n">kwarg1</span><span class="o">=</span><span class="n">val1</span><span class="p">,</span> <span class="n">kwarg2</span><span class="o">=</span><span class="n">val2</span><span class="p">,</span> <span class="o">...</span><span class="p">):</span>
    <span class="o">...</span>

<span class="n">solver</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">classname</span><span class="p">(</span><span class="n">f</span><span class="p">,</span>
                    <span class="n">f_args</span><span class="o">=</span><span class="p">[</span><span class="n">arg1</span><span class="p">,</span> <span class="n">arg2</span><span class="p">,</span> <span class="n">arg3</span><span class="p">,</span> <span class="o">...</span><span class="p">],</span>
                    <span class="n">f_kwargs</span><span class="o">=</span><span class="nb">dict</span><span class="p">(</span><span class="n">kwarg1</span><span class="o">=</span><span class="n">val1</span><span class="p">,</span> <span class="n">kwarg2</span><span class="o">=</span><span class="n">val2</span><span class="p">,</span> <span class="o">...</span><span class="p">))</span>

<span class="c"># Alternative setting of f_args and f_kwargs</span>
<span class="n">solver</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="n">f_args</span><span class="o">=</span><span class="p">[</span><span class="n">arg1</span><span class="p">,</span> <span class="n">arg2</span><span class="p">,</span> <span class="n">arg3</span><span class="p">,</span> <span class="o">...</span><span class="p">],</span>
           <span class="n">f_kwargs</span><span class="o">=</span><span class="nb">dict</span><span class="p">(</span><span class="n">kwarg1</span><span class="o">=</span><span class="n">val1</span><span class="p">,</span> <span class="n">kwarg2</span><span class="o">=</span><span class="n">val2</span><span class="p">,</span> <span class="o">...</span><span class="p">))</span>
</pre></div>
</div>
<p>Solvers will call <tt class="docutils literal"><span class="pre">f</span></tt> as <tt class="docutils literal"><span class="pre">f(u,</span> <span class="pre">t,</span> <span class="pre">*f_args,</span> <span class="pre">**f_kwargs)</span></tt>.</p>
</div>
<div class="section" id="termination-criterion-for-the-simulation">
<h3>Termination Criterion for the Simulation<a class="headerlink" href="#termination-criterion-for-the-simulation" title="Permalink to this headline">¶</a></h3>
<p>We know that the solution <span class="math">\(u\)</span> of the logistic equation approaches
<span class="math">\(R\)</span> as <span class="math">\(t\rightarrow\infty\)</span>. Instead of
using a trial and error process for determining an appropriate
time integral for integration, the <tt class="docutils literal"><span class="pre">solver.solve</span></tt> method accepts
a user-defined function <tt class="docutils literal"><span class="pre">terminate</span></tt> that can be used to implement
a criterion for terminating the solution process.
Mathematically, the relevant criterion is
<span class="math">\(||u-R||&lt;\hbox{tol}\)</span>, where tol is an acceptable
tolerance, say <span class="math">\(100\)</span> in the present case where <span class="math">\(R=10^5\)</span>.
The <tt class="docutils literal"><span class="pre">terminate</span></tt> function implements the criterion
and returns true if the criterion is met:</p>
<div class="highlight-python"><pre>def terminate(u, t, step_no):
    """u[step_no] holds (the most recent) solution at t[step_no]."""
    return abs(u[step_no] - R) &lt; tol:</pre>
</div>
<p>Note that the simulation is anyway stopped for <span class="math">\(t &gt; T\)</span> so <span class="math">\(T\)</span>
must be large enough for the termination criterion to be reached (if not,
a warning will be issued).
With a <tt class="docutils literal"><span class="pre">terminate</span></tt> function it is also convenient to specify the
time step <tt class="docutils literal"><span class="pre">dt</span></tt> and not the total number of time steps.</p>
<p>A complete program can be as follows (<a class="reference external" href="https://github.com/hplgit/odespy/blob/master/doc/src/odespy/src-odespy/logistic5.py">logistic5.py</a>):</p>
<div class="highlight-python"><pre>def f(u, t):
    return a*u*(1 - u/R)

a = 2
R = 1E+5
A = 1

import odespy, numpy
solver = odespy.RK4(f)
solver.set_initial_condition(A)

T = 20  # end of simulation
dt = 0.25
N = int(round(T/dt))
time_points = numpy.linspace(0, T, N+1)

tol = 100   # tolerance for termination criterion

def terminate(u, t, step_no):
    """u[step_no] holds (the most recent) solution at t[step_no]."""
    return abs(u[step_no] - R) &lt; tol:

u, t = solver.solve(time_points, terminate)
print 'Final u(t=%g)=%g' % (t[-1], u[-1])

from matplotlib.pyplot import *
plot(t, u, 'r-')
savefig('tmp.png')
show()</pre>
</div>
</div>
<div class="section" id="a-class-based-implementation">
<h3>A Class-Based Implementation<a class="headerlink" href="#a-class-based-implementation" title="Permalink to this headline">¶</a></h3>
<p>The previous code example can be recast into a more class-based
(&#8220;object-oriented programming&#8221;) example. We lump all data related
to the problem (the &#8220;physics&#8221;) into a problem class <tt class="docutils literal"><span class="pre">Logistic</span></tt>, while
all data related to the numerical solution and its quality are
taken care of by class <tt class="docutils literal"><span class="pre">Solver</span></tt>. The code below illustrates
the ideas (<a class="reference external" href="https://github.com/hplgit/odespy/blob/master/doc/src/odespy/src-odespy/logistic26py">logistic6.py</a>):</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">np</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="kn">as</span> <span class="nn">mpl</span>
<span class="kn">import</span> <span class="nn">odespy</span>

<span class="k">class</span> <span class="nc">Logistic</span><span class="p">:</span>
    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">R</span><span class="p">,</span> <span class="n">A</span><span class="p">,</span> <span class="n">T</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        `a` is (initial growth rate), `R` the carrying capacity,</span>
<span class="sd">        `A` the initial amount of u, and `T` is some (very) total</span>
<span class="sd">        simulation time when `u` is very close to the asymptotic</span>
<span class="sd">        value `R`.</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">R</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">A</span> <span class="o">=</span> <span class="n">a</span><span class="p">,</span> <span class="n">R</span><span class="p">,</span> <span class="n">A</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">tol</span> <span class="o">=</span> <span class="mf">0.01</span><span class="o">*</span><span class="n">R</span> <span class="c"># tolerance for termination criterion</span>

    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;Right-hand side of the ODE.&quot;&quot;&quot;</span>
        <span class="n">a</span><span class="p">,</span> <span class="n">R</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">R</span>  <span class="c"># short form</span>
        <span class="k">return</span> <span class="n">a</span><span class="o">*</span><span class="n">u</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u</span><span class="o">/</span><span class="n">R</span><span class="p">)</span>

    <span class="k">def</span> <span class="nf">terminate</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">step_no</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;u[step_no] holds solution at t[step_no].&quot;&quot;&quot;</span>
        <span class="k">return</span> <span class="nb">abs</span><span class="p">(</span><span class="n">u</span><span class="p">[</span><span class="n">step_no</span><span class="p">]</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">R</span><span class="p">)</span> <span class="o">&lt;</span> <span class="bp">self</span><span class="o">.</span><span class="n">tol</span>

    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">a</span><span class="p">,</span> <span class="n">R</span><span class="p">,</span> <span class="n">A</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">R</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">A</span>  <span class="c"># short form</span>
        <span class="k">return</span> <span class="n">R</span><span class="o">*</span><span class="n">A</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="o">*</span><span class="n">t</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">R</span> <span class="o">+</span> <span class="n">A</span><span class="o">*</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="o">*</span><span class="n">t</span><span class="p">)</span> <span class="o">-</span> <span class="mi">1</span><span class="p">))</span>


<span class="k">class</span> <span class="nc">Solver</span><span class="p">:</span>
    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">problem</span><span class="p">,</span> <span class="n">dt</span><span class="p">,</span> <span class="n">method</span><span class="o">=</span><span class="s">&#39;RK4&#39;</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">problem</span> <span class="o">=</span> <span class="n">problem</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">dt</span> <span class="o">=</span> <span class="n">dt</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">method_class</span> <span class="o">=</span> <span class="nb">eval</span><span class="p">(</span><span class="s">&#39;odespy.&#39;</span> <span class="o">+</span> <span class="n">method</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">N</span> <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="nb">round</span><span class="p">(</span><span class="n">T</span><span class="o">/</span><span class="n">dt</span><span class="p">))</span>

    <span class="k">def</span> <span class="nf">solve</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">solver</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">method_class</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">problem</span><span class="o">.</span><span class="n">f</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">solver</span><span class="o">.</span><span class="n">set_initial_condition</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">problem</span><span class="o">.</span><span class="n">A</span><span class="p">)</span>
        <span class="n">time_points</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">problem</span><span class="o">.</span><span class="n">T</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">N</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">u</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span>
            <span class="n">time_points</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">problem</span><span class="o">.</span><span class="n">terminate</span><span class="p">)</span>
        <span class="k">print</span> <span class="s">&#39;Final u(t=</span><span class="si">%g</span><span class="s">)=</span><span class="si">%g</span><span class="s">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">t</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">],</span> <span class="n">u</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">])</span>

    <span class="k">def</span> <span class="nf">plot</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="n">mpl</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">u</span><span class="p">,</span> <span class="s">&#39;r-&#39;</span><span class="p">,</span>
                 <span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="p">),</span> <span class="s">&#39;bo&#39;</span><span class="p">)</span>
        <span class="n">mpl</span><span class="o">.</span><span class="n">legend</span><span class="p">([</span><span class="s">&#39;numerical&#39;</span><span class="p">,</span> <span class="s">&#39;exact&#39;</span><span class="p">])</span>
        <span class="n">mpl</span><span class="o">.</span><span class="n">savefig</span><span class="p">(</span><span class="s">&#39;tmp.png&#39;</span><span class="p">)</span>
        <span class="n">mpl</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>

<span class="k">def</span> <span class="nf">main</span><span class="p">():</span>
    <span class="n">problem</span> <span class="o">=</span> <span class="n">Logistic</span><span class="p">(</span><span class="n">a</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span> <span class="n">R</span><span class="o">=</span><span class="mf">1E+5</span><span class="p">,</span> <span class="n">A</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">T</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
    <span class="n">solver</span> <span class="o">=</span> <span class="n">Solver</span><span class="p">(</span><span class="n">problem</span><span class="p">,</span> <span class="n">dt</span><span class="o">=</span><span class="mf">0.25</span><span class="p">,</span> <span class="n">method</span><span class="o">=</span><span class="s">&#39;RK4&#39;</span><span class="p">)</span>
    <span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">()</span>
    <span class="n">solver</span><span class="o">.</span><span class="n">plot</span><span class="p">()</span>

<span class="k">if</span> <span class="n">__name__</span> <span class="o">==</span> <span class="s">&#39;__main__&#39;</span><span class="p">:</span>
    <span class="n">main</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="section" id="using-other-symbols">
<h3>Using Other Symbols<a class="headerlink" href="#using-other-symbols" title="Permalink to this headline">¶</a></h3>
<p>The Odespy package applies <tt class="docutils literal"><span class="pre">u</span></tt> for the unknown function or vector of
unknown functions and <tt class="docutils literal"><span class="pre">t</span></tt> as the name of the independent
variable. Many problems involve other symbols for functions and
independent variables. These symbols should be reflected in the user&#8217;s
code.  For example, here is a coding example involving the logistic
equation written as <span class="math">\(y'(x)=au(x)(1-u(x)/R(x))\)</span>, where now a variable
<span class="math">\(R=R(x)\)</span> is considered. Following the setup from the very first
program above solving the logistic ODE, we can easily introduce our
own nomenclature (<a class="reference external" href="https://github.com/hplgit/odespy/blob/master/doc/src/odespy/src-odespy/logistic7.py">logistic7.py</a>):</p>
<div class="highlight-python"><div class="highlight"><pre><span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="n">y</span><span class="p">,</span> <span class="n">x</span><span class="p">):</span>
    <span class="k">return</span> <span class="n">a</span><span class="o">*</span><span class="n">y</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">y</span><span class="o">/</span><span class="n">R</span><span class="p">)</span>

<span class="n">a</span> <span class="o">=</span> <span class="mi">2</span><span class="p">;</span>  <span class="n">R</span> <span class="o">=</span> <span class="mf">1E+5</span><span class="p">;</span>  <span class="n">A</span> <span class="o">=</span> <span class="mi">1</span>

<span class="kn">import</span> <span class="nn">odespy</span><span class="o">,</span> <span class="nn">numpy</span>
<span class="n">solver</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">RK4</span><span class="p">(</span><span class="n">f</span><span class="p">)</span>
<span class="n">solver</span><span class="o">.</span><span class="n">set_initial_condition</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>

<span class="n">L</span> <span class="o">=</span> <span class="mi">10</span>  <span class="c"># end of x domain</span>
<span class="n">N</span> <span class="o">=</span> <span class="mi">30</span>  <span class="c"># no of time steps</span>
<span class="n">x_points</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">L</span><span class="p">,</span> <span class="n">N</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>
<span class="n">y</span><span class="p">,</span> <span class="n">x</span> <span class="o">=</span> <span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">x_points</span><span class="p">)</span>

<span class="kn">from</span> <span class="nn">matplotlib.pyplot</span> <span class="kn">import</span> <span class="o">*</span>
<span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="s">&#39;r-&#39;</span><span class="p">)</span>
<span class="n">xlabel</span><span class="p">(</span><span class="s">&#39;x&#39;</span><span class="p">);</span> <span class="n">ylabel</span><span class="p">(</span><span class="s">&#39;y&#39;</span><span class="p">)</span>
<span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
<p>As shown, we use <tt class="docutils literal"><span class="pre">y</span></tt> for <tt class="docutils literal"><span class="pre">u</span></tt>, <tt class="docutils literal"><span class="pre">x</span></tt> for <tt class="docutils literal"><span class="pre">t</span></tt>, and <tt class="docutils literal"><span class="pre">x_points</span></tt> instead
of <tt class="docutils literal"><span class="pre">time_points</span></tt>.</p>
</div>
<div class="section" id="example-solving-an-ode-system">
<h3>Example: Solving an ODE System<a class="headerlink" href="#example-solving-an-ode-system" title="Permalink to this headline">¶</a></h3>
<p>We shall now explain how to solve a system of ODEs using a scalar
second-order ODE as starting point.
The angle <span class="math">\(\theta\)</span> of a pendulum with mass <span class="math">\(m\)</span> and length <span class="math">\(L\)</span>
is governed by the equation
(neglecting air resistance for simplicity)</p>
<div class="math">
\[mL\ddot\theta + mg\sin\theta = 0,\quad \theta (0)=\Theta,\
\dot\theta (0)=0 .\]</div>
<p>A dot over <span class="math">\(\theta\)</span> implies differentiation with respect to time.
The ODE can be written as <span class="math">\(\ddot\theta + c\sin\theta=0\)</span> by
introducing <span class="math">\(c = g/L\)</span>.
This problem must be expressed as a first-order ODE system if it is
going to be solved by the tools in the Odespy package.
Introducing <span class="math">\(\omega = \dot\theta\)</span> (the angular velocity) as auxiliary
unknown, we get the system</p>
<div class="math">
\[\begin{split}\dot\theta &amp;= \omega,\\
\dot\omega &amp;= -c\sin\theta,\end{split}\]</div>
<p>with <span class="math">\(\theta(0)=\Theta\)</span> and <span class="math">\(\omega(0)=0\)</span>.</p>
<p>Now the <tt class="docutils literal"><span class="pre">f</span></tt> function must return a list or array with the two
right-hand side functions:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
    <span class="n">theta</span><span class="p">,</span> <span class="n">omega</span> <span class="o">=</span> <span class="n">u</span>
    <span class="k">return</span> <span class="p">[</span><span class="n">omega</span><span class="p">,</span> <span class="o">-</span><span class="n">c</span><span class="o">*</span><span class="n">sin</span><span class="p">(</span><span class="n">theta</span><span class="p">)]</span>
</pre></div>
</div>
<p>Note that when we have a system of ODEs with <tt class="docutils literal"><span class="pre">n</span></tt> components, the <tt class="docutils literal"><span class="pre">u</span></tt>
object sent to the <tt class="docutils literal"><span class="pre">f</span></tt> function is an array of length <tt class="docutils literal"><span class="pre">n</span></tt>,
representing the value of all components in the ODE system at time <tt class="docutils literal"><span class="pre">t</span></tt>.
Here we extract the two components of <tt class="docutils literal"><span class="pre">u</span></tt> in separate local variables
with names equal to what is used in the mathematical description of
the current problem.</p>
<p>The initial conditions must be specified as a list:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">solver</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">Heun</span><span class="p">(</span><span class="n">f</span><span class="p">)</span>
<span class="n">solver</span><span class="o">.</span><span class="n">set_initial_condition</span><span class="p">([</span><span class="n">Theta</span><span class="p">,</span> <span class="mi">0</span><span class="p">])</span>
</pre></div>
</div>
<p>To specify the time points we here first decide on a number of periods
(oscillations back and forth) to simulate and then on the time resolution
of each period. Note that when <span class="math">\(\Theta\)</span> is small we can replace
<span class="math">\(\sin\theta\)</span> by <span class="math">\(\theta\)</span> and find an analytical
solution
<span class="math">\(\theta (t)=\Theta\cos\left(\sqrt{c}t\right)\)</span>
whose period is <span class="math">\(2\pi/\sqrt{c}\)</span>. We use this expression
as an approximation for the period also when <span class="math">\(\Theta\)</span> is not
small.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">freq</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">c</span><span class="p">)</span>      <span class="c"># frequency of oscillations when Theta is small</span>
<span class="n">period</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="o">/</span><span class="n">freq</span>  <span class="c"># the period of the oscillations</span>
<span class="n">T</span> <span class="o">=</span> <span class="mi">10</span><span class="o">*</span><span class="n">period</span>       <span class="c"># final time</span>
<span class="n">N_per_period</span> <span class="o">=</span> <span class="mi">20</span>   <span class="c"># resolution of one period</span>
<span class="n">N</span> <span class="o">=</span> <span class="n">N_per_period</span><span class="o">*</span><span class="n">period</span>
<span class="n">time_points</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">T</span><span class="p">,</span> <span class="n">N</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>

<span class="n">u</span><span class="p">,</span> <span class="n">t</span> <span class="o">=</span> <span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">time_points</span><span class="p">)</span>
</pre></div>
</div>
<p>The <tt class="docutils literal"><span class="pre">u</span></tt> returned from <tt class="docutils literal"><span class="pre">solver.solve</span></tt> is a two-dimensional array, where the
columns hold the various solution functions of the ODE system. That is,
the first column holds <span class="math">\(\theta\)</span> and the second column holds
<span class="math">\(\omega\)</span>. For convenience we extract the individual solution
components in individual arrays:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">theta</span> <span class="o">=</span> <span class="n">u</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span>
<span class="n">omega</span> <span class="o">=</span> <span class="n">u</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]</span>

<span class="kn">from</span> <span class="nn">matplotlib.pyplot</span> <span class="kn">import</span> <span class="o">*</span>
<span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span> <span class="n">theta</span><span class="p">,</span> <span class="s">&#39;r-&#39;</span><span class="p">)</span>
<span class="n">savefig</span><span class="p">(</span><span class="s">&#39;tmp.png&#39;</span><span class="p">)</span>
<span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
<p>The complete program is available in the file <a class="reference external" href="https://github.com/hplgit/odespy/blob/master/doc/src/odespy/src-odespy/osc1a.py">osc1a.py</a>.</p>
<p>Looking at the plot reveals that the numerical solution has
an alarming feature: the amplitude grows (indicating increasing
energy in the system). Changing <tt class="docutils literal"><span class="pre">T</span></tt> to 28 periods instead of 10
makes the numerical solution explode.
The increasing amplitude is a numerical artifact that some of the simple
solution methods suffer from.</p>
<div class="figure">
<img alt="_images/osc1a.png" src="_images/osc1a.png" style="width: 500px;" />
<p class="caption">Heun&#8217;s method used to simulate oscillations of a pendulum</p>
</div>
<p>Using a more sophisticated method, say the 4-th order Runge-Kutta method,
is just a matter of substituting <tt class="docutils literal"><span class="pre">Heun</span></tt> by <tt class="docutils literal"><span class="pre">RK4</span></tt>:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">solver</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">RK4</span><span class="p">(</span><span class="n">f</span><span class="p">)</span>
<span class="n">solver</span><span class="o">.</span><span class="n">set_initial_condition</span><span class="p">([</span><span class="n">Theta</span><span class="p">,</span> <span class="mi">0</span><span class="p">])</span>
<span class="n">freq</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">c</span><span class="p">)</span>      <span class="c"># frequency of oscillations when Theta is small</span>
<span class="n">period</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="o">/</span><span class="n">freq</span>  <span class="c"># the period of the oscillations</span>
<span class="n">T</span> <span class="o">=</span> <span class="mi">10</span><span class="o">*</span><span class="n">period</span>       <span class="c"># final time</span>
<span class="n">N_per_period</span> <span class="o">=</span> <span class="mi">20</span>   <span class="c"># resolution of one period</span>
<span class="n">N</span> <span class="o">=</span> <span class="n">N_per_period</span><span class="o">*</span><span class="n">period</span>
<span class="n">time_points</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">T</span><span class="p">,</span> <span class="n">N</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>

<span class="n">u</span><span class="p">,</span> <span class="n">t</span> <span class="o">=</span> <span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">time_points</span><span class="p">)</span>

<span class="n">theta</span> <span class="o">=</span> <span class="n">u</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span>
<span class="n">omega</span> <span class="o">=</span> <span class="n">u</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]</span>

<span class="kn">from</span> <span class="nn">matplotlib.pyplot</span> <span class="kn">import</span> <span class="o">*</span>
<span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span> <span class="n">theta</span><span class="p">,</span> <span class="s">&#39;r-&#39;</span><span class="p">)</span>
<span class="n">savefig</span><span class="p">(</span><span class="s">&#39;tmp.png&#39;</span><span class="p">)</span>
<span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
<p>The amplitude now becomes constant in time as expected.</p>
</div>
<div class="section" id="testing-several-methods">
<h3>Testing Several Methods<a class="headerlink" href="#testing-several-methods" title="Permalink to this headline">¶</a></h3>
<p>We shall now make a more advanced solver by
extending the previous example. More specifically, we shall</p>
<blockquote>
<div><ul class="simple">
<li>represent the right-hand side function as class,</li>
<li>compare several different solvers,</li>
<li>compute error of numerical solutions.</li>
</ul>
</div></blockquote>
<p>Since we want to compare numerical errors in the various
solvers we need a test problem where the exact solution is known.
Approximating <span class="math">\(\sin(\theta)\)</span> by <span class="math">\(\theta\)</span>
(valid for small <span class="math">\(\theta\)</span>), gives the ODE system</p>
<div class="math">
\[\begin{split}\dot\theta &amp;= \omega,\\
\dot\omega &amp;= -c\theta,\end{split}\]</div>
<p>with <span class="math">\(\theta(0)=\Theta\)</span> and <span class="math">\(\omega(0)=0\)</span>.</p>
<p>Right-hand side functions with parameters can be handled by
including extra arguments via the <tt class="docutils literal"><span class="pre">f_args</span></tt> and <tt class="docutils literal"><span class="pre">f_kwargs</span></tt> functionality,
or by using a class where the parameters are attributes and
an <tt class="docutils literal"><span class="pre">f</span></tt> method defines <span class="math">\(f(u,t)\)</span>. The section <a class="reference internal" href="#ode-sec-exgr-farg"><em>Parameters in the Right-Hand Side Function</em></a> exemplifies the details.
A minimal class representation of the right-hand side
function in the present case looks like this:</p>
<div class="highlight-py"><div class="highlight"><pre><span class="k">class</span> <span class="nc">Problem</span><span class="p">:</span>
    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">c</span><span class="p">,</span> <span class="n">Theta</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">Theta</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">c</span><span class="p">),</span> <span class="nb">float</span><span class="p">(</span><span class="n">Theta</span><span class="p">)</span>

    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">theta</span><span class="p">,</span> <span class="n">omega</span> <span class="o">=</span> <span class="n">u</span><span class="p">;</span>  <span class="n">c</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">c</span>
        <span class="k">return</span> <span class="p">[</span><span class="n">omega</span><span class="p">,</span> <span class="o">-</span><span class="n">c</span><span class="o">*</span><span class="n">theta</span><span class="p">]</span>

<span class="n">problem</span> <span class="o">=</span> <span class="n">Problem</span><span class="p">(</span><span class="n">c</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">Theta</span><span class="o">=</span><span class="n">pi</span><span class="o">/</span><span class="mi">4</span><span class="p">)</span>
</pre></div>
</div>
<p>It would be convenient to add an attribute <tt class="docutils literal"><span class="pre">period</span></tt> which holds
an estimate of the period of oscillations as we need this for
deciding on the complete time interval for solving the differential
equations. An appropriate extension of class <tt class="docutils literal"><span class="pre">Problem</span></tt> is therefore</p>
<div class="highlight-python"><div class="highlight"><pre><span class="k">class</span> <span class="nc">Problem</span><span class="p">:</span>
    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">c</span><span class="p">,</span> <span class="n">Theta</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">Theta</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">c</span><span class="p">),</span> <span class="nb">float</span><span class="p">(</span><span class="n">Theta</span><span class="p">)</span>

        <span class="bp">self</span><span class="o">.</span><span class="n">freq</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">c</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">period</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">freq</span>

    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">theta</span><span class="p">,</span> <span class="n">omega</span> <span class="o">=</span> <span class="n">u</span><span class="p">;</span>  <span class="n">c</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">c</span>
        <span class="k">return</span> <span class="p">[</span><span class="n">omega</span><span class="p">,</span> <span class="o">-</span><span class="n">c</span><span class="o">*</span><span class="n">theta</span><span class="p">]</span>

<span class="n">problem</span> <span class="o">=</span> <span class="n">Problem</span><span class="p">(</span><span class="n">c</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">Theta</span><span class="o">=</span><span class="n">pi</span><span class="o">/</span><span class="mi">4</span><span class="p">)</span>
</pre></div>
</div>
<p>The second extension is to loop over many solvers. All
solvers can be listed by</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="kn">import</span> <span class="nn">odespy</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">methods</span> <span class="o">=</span> <span class="n">list_all_solvers</span><span class="p">()</span>
<span class="gp">&gt;&gt;&gt; </span><span class="k">for</span> <span class="n">method</span> <span class="ow">in</span> <span class="n">methods</span><span class="p">:</span>
<span class="gp">... </span>  <span class="k">print</span> <span class="n">method</span>
<span class="gp">...</span>
<span class="go">AdamsBashMoulton2</span>
<span class="go">AdamsBashMoulton3</span>
<span class="go">AdamsBashforth2</span>
<span class="gp">...</span>
<span class="go">Vode</span>
<span class="go">lsoda_scipy</span>
<span class="go">odefun_sympy</span>
<span class="go">odelab</span>
</pre></div>
</div>
<p>A similar function, <tt class="docutils literal"><span class="pre">list_available_solvers</span></tt>, returns a list of the
names of the solvers that are available in the current installation
(e.g., the <tt class="docutils literal"><span class="pre">Vode</span></tt> solver is only available if the comprehensive
<tt class="docutils literal"><span class="pre">scipy</span></tt> package is installed).
This is the list that is usually most relevant.</p>
<p>For now we explicitly choose a subset of the commonly available solvers:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">odespy</span>
<span class="n">solvers</span> <span class="o">=</span> <span class="p">[</span>
    <span class="n">odespy</span><span class="o">.</span><span class="n">ThetaRule</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">f</span><span class="p">,</span> <span class="n">theta</span><span class="o">=</span><span class="mi">0</span><span class="p">),</span>   <span class="c"># Forward Euler</span>
    <span class="n">odespy</span><span class="o">.</span><span class="n">ThetaRule</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">f</span><span class="p">,</span> <span class="n">theta</span><span class="o">=</span><span class="mf">0.5</span><span class="p">),</span> <span class="c"># Midpoint method</span>
    <span class="n">odespy</span><span class="o">.</span><span class="n">ThetaRule</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">f</span><span class="p">,</span> <span class="n">theta</span><span class="o">=</span><span class="mi">1</span><span class="p">),</span>   <span class="c"># Backward Euler</span>
    <span class="n">odespy</span><span class="o">.</span><span class="n">RK4</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">f</span><span class="p">),</span>
    <span class="n">odespy</span><span class="o">.</span><span class="n">MidpointIter</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">f</span><span class="p">,</span> <span class="n">max_iter</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span> <span class="n">eps_iter</span><span class="o">=</span><span class="mf">0.01</span><span class="p">),</span>
    <span class="n">odespy</span><span class="o">.</span><span class="n">LeapfrogFiltered</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">f</span><span class="p">),</span>
    <span class="p">]</span>
</pre></div>
</div>
<p>To see what a method is and its arguments to the constructor, invoke
the doc string of the class, e.g., <tt class="docutils literal"><span class="pre">help(ThetaRule)</span></tt> inside a
Python shell like IPython, or run <tt class="docutils literal"><span class="pre">pydoc</span> <span class="pre">odespy.ThetaRule</span></tt> in a
terminal window, or invoke the <a class="reference external" href="http://hplgit.github.com/odespy/doc/api/">Odespy API documentation</a>.</p>
<p>It will be evident that the <tt class="docutils literal"><span class="pre">ThetaRule</span></tt> solver with <tt class="docutils literal"><span class="pre">theta=0</span></tt> and
<tt class="docutils literal"><span class="pre">theta=1</span></tt> (Forward and Backward Euler methods) gives growing and
decaying amplitudes, respectively, while the other solvers are
capable of reproducing the constant amplitude of the oscillations of
in the current mathematical model.</p>
<p>The loop over the chosen solvers may look like</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">N_per_period</span> <span class="o">=</span> <span class="mi">20</span>
<span class="n">T</span> <span class="o">=</span> <span class="mi">3</span><span class="o">*</span><span class="n">problem</span><span class="o">.</span><span class="n">period</span>   <span class="c"># final time</span>
<span class="kn">import</span> <span class="nn">numpy</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="kn">as</span> <span class="nn">mpl</span>
<span class="n">legends</span> <span class="o">=</span> <span class="p">[]</span>

<span class="k">for</span> <span class="n">solver</span> <span class="ow">in</span> <span class="n">solvers</span><span class="p">:</span>
    <span class="n">solver_name</span> <span class="o">=</span> <span class="nb">str</span><span class="p">(</span><span class="n">solver</span><span class="p">)</span>  <span class="c"># short description of solver</span>
    <span class="k">print</span> <span class="n">solver_name</span>

    <span class="n">solver</span><span class="o">.</span><span class="n">set_initial_condition</span><span class="p">([</span><span class="n">problem</span><span class="o">.</span><span class="n">Theta</span><span class="p">,</span> <span class="mi">0</span><span class="p">])</span>
    <span class="n">N</span> <span class="o">=</span> <span class="n">N_per_period</span><span class="o">*</span><span class="n">problem</span><span class="o">.</span><span class="n">period</span>
    <span class="n">time_points</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">T</span><span class="p">,</span> <span class="n">N</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>

    <span class="n">u</span><span class="p">,</span> <span class="n">t</span> <span class="o">=</span> <span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">time_points</span><span class="p">)</span>

    <span class="n">theta</span> <span class="o">=</span> <span class="n">u</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">legends</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">solver_name</span><span class="p">)</span>
    <span class="n">mpl</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">t</span><span class="p">,</span> <span class="n">theta</span><span class="p">)</span>
    <span class="n">mpl</span><span class="o">.</span><span class="n">hold</span><span class="p">(</span><span class="s">&#39;on&#39;</span><span class="p">)</span>
<span class="n">mpl</span><span class="o">.</span><span class="n">legend</span><span class="p">(</span><span class="n">legends</span><span class="p">)</span>
<span class="n">mpl</span><span class="o">.</span><span class="n">savefig</span><span class="p">(</span><span class="n">__file__</span><span class="p">[:</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span> <span class="o">+</span> <span class="s">&#39;.png&#39;</span><span class="p">)</span>
<span class="n">mpl</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
<p>A complete program is available as <a class="reference external" href="https://github.com/hplgit/odespy/blob/master/doc/src/odespy/src-odespy/osc2.py">osc2.py</a>.</p>
<div class="figure">
<img alt="_images/osc2.png" src="_images/osc2.png" style="width: 500px;" />
<p class="caption">Comparison of methods for solving the ODE system for a pendulum</p>
</div>
</div>
<div class="section" id="make-a-subclass-of-problem">
<h3>Make a Subclass of Problem<a class="headerlink" href="#make-a-subclass-of-problem" title="Permalink to this headline">¶</a></h3>
<p>Odespy features a module <tt class="docutils literal"><span class="pre">problems</span></tt> for defining ODE problems.
There is a superclass <tt class="docutils literal"><span class="pre">Problem</span></tt> in this module defining what we
expect of information about an ODE problem, as well as some
convenience functions that are inherited in subclasses.
A rough sketch of class <tt class="docutils literal"><span class="pre">Problem</span></tt> is listed here:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="k">class</span> <span class="nc">Problem</span><span class="p">:</span>
    <span class="n">stiff</span> <span class="o">=</span> <span class="bp">False</span>    <span class="c"># classification of the problem is stiff or not</span>
    <span class="n">complex_</span> <span class="o">=</span> <span class="bp">False</span> <span class="c"># True if f(u,t) is complex valued</span>
    <span class="n">not_suitable_solvers</span> <span class="o">=</span> <span class="p">[]</span>  <span class="c"># list solvers that should be be used</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&#39;&#39;</span>     <span class="c"># one-line problem description</span>

    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="k">pass</span>

    <span class="k">def</span> <span class="nf">__contains__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">attr</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;Return True if attr is a method in instance self.&quot;&quot;&quot;</span>

    <span class="k">def</span> <span class="nf">terminate</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">step_number</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;Default terminate function, always returning False.&quot;&quot;&quot;</span>
        <span class="k">return</span> <span class="bp">False</span>

    <span class="k">def</span> <span class="nf">default_parameters</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Compute suitable time_points, atol/rtol, etc. for the</span>
<span class="sd">        particular problem. Useful for quick generation of test</span>
<span class="sd">        cases, demos, unit tests, etc.</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="k">return</span> <span class="p">{}</span>

    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;Implementation of the exact solution.&quot;&quot;&quot;</span>
        <span class="k">return</span> <span class="bp">None</span>
</pre></div>
</div>
<p>Subclasses of <tt class="docutils literal"><span class="pre">Problem</span></tt> typically implements the constructor, for
registering parameters in the ODE and the initial condition, and
a method <tt class="docutils literal"><span class="pre">f</span></tt> for defining the right-hand side. For implicit solution method
we may provide a method <tt class="docutils literal"><span class="pre">jac</span></tt> returning the Jacobian of <span class="math">\(f(u,t)\)</span> with
respect to <span class="math">\(u\)</span>. Some problems may
also register an analytical solution in <tt class="docutils literal"><span class="pre">u_exact</span></tt>. Here is an
example of implementing the logistic ODE from the section <a class="reference internal" href="#ode-sec-exgr"><em>First Example: Logistic Growth</em></a>:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">odespy</span>

<span class="k">class</span> <span class="nc">Logistic</span><span class="p">(</span><span class="n">odespy</span><span class="o">.</span><span class="n">problems</span><span class="o">.</span><span class="n">Problem</span><span class="p">):</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&quot;Logistic equation&quot;</span>

    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">R</span><span class="p">,</span> <span class="n">A</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">a</span> <span class="o">=</span> <span class="n">a</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">R</span> <span class="o">=</span> <span class="n">R</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">U0</span> <span class="o">=</span> <span class="n">A</span>

    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">a</span><span class="p">,</span> <span class="n">R</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">R</span>  <span class="c"># short form</span>
        <span class="k">return</span> <span class="n">a</span><span class="o">*</span><span class="n">u</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u</span><span class="o">/</span><span class="n">R</span><span class="p">)</span>

    <span class="k">def</span> <span class="nf">jac</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">a</span><span class="p">,</span> <span class="n">R</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">R</span>  <span class="c"># short form</span>
        <span class="k">return</span> <span class="n">a</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u</span><span class="o">/</span><span class="n">R</span><span class="p">)</span> <span class="o">+</span> <span class="n">a</span><span class="o">*</span><span class="n">u</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="mf">1.</span><span class="o">/</span><span class="n">R</span><span class="p">)</span>

    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">a</span><span class="p">,</span> <span class="n">R</span><span class="p">,</span> <span class="n">U0</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">R</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">U0</span>  <span class="c"># short form</span>
        <span class="k">return</span> <span class="n">R</span><span class="o">*</span><span class="n">U0</span><span class="o">*</span><span class="n">numpy</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="o">*</span><span class="n">t</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">R</span> <span class="o">+</span> <span class="n">U0</span><span class="o">*</span><span class="p">(</span><span class="n">numpy</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="o">*</span><span class="n">t</span><span class="p">)</span> <span class="o">-</span> <span class="mi">1</span><span class="p">))</span>
</pre></div>
</div>
<p>The <tt class="docutils literal"><span class="pre">stiff</span></tt>, <tt class="docutils literal"><span class="pre">complex_</span></tt>, and <tt class="docutils literal"><span class="pre">not_suitable_solvers</span></tt> class variables
can just be inherited. Note that <tt class="docutils literal"><span class="pre">u_exact</span></tt> should work for a vector <tt class="docutils literal"><span class="pre">t</span></tt>
so <tt class="docutils literal"><span class="pre">numpy</span></tt> versions of mathematical functions must be used.</p>
<p>The initial condition is by convention stored as the attribute <tt class="docutils literal"><span class="pre">U0</span></tt>
in a subclass of <tt class="docutils literal"><span class="pre">Problem</span></tt>, and specified as argument to the constructor.</p>
<p>Here are the typical steps when using such a problem class:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">problem</span> <span class="o">=</span> <span class="n">Logistic</span><span class="p">(</span><span class="n">a</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span> <span class="n">R</span><span class="o">=</span><span class="mf">1E+5</span><span class="p">,</span> <span class="n">A</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
<span class="n">solver</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">RK4</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">f</span><span class="p">)</span>
<span class="n">solver</span><span class="o">.</span><span class="n">set_initial_condition</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">U0</span><span class="p">)</span>
<span class="n">u</span><span class="p">,</span> <span class="n">t</span> <span class="o">=</span> <span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">time_points</span><span class="p">)</span>
</pre></div>
</div>
<p>The problem class may also feature additional methods:</p>
<div class="highlight-python"><pre>class MyProblem(odespy.problems.Problem)
    ...
    def constraints(self, u, t):
        """Python function for additional constraints: g(u,t)=0."""

    def define_command_line_arguments(self, parser):
        """
        Initialize an argparse object for reading command-line
        option-value pairs. `parser` is an ``argparse`` object.
        """

    def verify(self, u, t, atol=None, rtol=None):
        """
        Return True if u at time points t coincides with an exact
        solution within the prescribed tolerances. If one of the
        tolerances is None, return max computed error (infinity
        norm). Return None if the solution cannot be verified.
        """</pre>
</div>
<p>The module <tt class="docutils literal"><span class="pre">odespy.problems</span></tt> contains many predefined ODE problems.</p>
</div>
<div class="section" id="using-adaptive-methods">
<h3>Using Adaptive Methods<a class="headerlink" href="#using-adaptive-methods" title="Permalink to this headline">¶</a></h3>
<p>The solvers used in the previous examples have all employed a constant
time step <span class="math">\(\Delta t\)</span>. Many solvers available through the Odespy
interface are adaptive in the sense that <span class="math">\(\Delta t\)</span> is adjusted
throughout the solution process to meet a prescribed tolerance for
the estimated error.</p>
<p>Simple methods such as <tt class="docutils literal"><span class="pre">RK4</span></tt> apply time steps</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">dt</span> <span class="o">=</span> <span class="n">time_points</span><span class="p">[</span><span class="n">k</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">time_points</span><span class="p">[</span><span class="n">k</span><span class="p">]</span>
</pre></div>
</div>
<p>while adaptive methods will use several (smaller) time steps than <tt class="docutils literal"><span class="pre">dt</span></tt>
in each <tt class="docutils literal"><span class="pre">dt</span></tt> interval to ensure that the estimated numerical error is
smaller than some prescribed tolerance. The estimated numerical error
may be a rather crude quantitative measure of the true numerical
error (which we do not know since the exact solution of the problem is
in general not known).</p>
<p>Some adaptive solvers record the intermediate solutions in each <tt class="docutils literal"><span class="pre">dt</span></tt>
interval in arrays <tt class="docutils literal"><span class="pre">self.u_all</span></tt> and <tt class="docutils literal"><span class="pre">self.t_all</span></tt>.  Examples include
<tt class="docutils literal"><span class="pre">RKFehlberg</span></tt>, <tt class="docutils literal"><span class="pre">Fehlberg</span></tt>, <tt class="docutils literal"><span class="pre">DormandPrince</span></tt>, <tt class="docutils literal"><span class="pre">CashKarp</span></tt>, and
<tt class="docutils literal"><span class="pre">BogackiShampine</span></tt>. Other adaptive solvers (<tt class="docutils literal"><span class="pre">Vode</span></tt>, <tt class="docutils literal"><span class="pre">Lsode</span></tt>, <tt class="docutils literal"><span class="pre">Lsoda</span></tt>,
<tt class="docutils literal"><span class="pre">RKC</span></tt>, <tt class="docutils literal"><span class="pre">RKF45</span></tt>, etc.)  do not give access to intermediate solution
steps between the user-given time points, specified in the
<tt class="docutils literal"><span class="pre">solver.solve</span></tt> call, and then we only have access to the solution at
the user-given time points as returned by this call. One can run <tt class="docutils literal"><span class="pre">if</span>
<span class="pre">solver.has_u_t_all()</span></tt> to test if the <tt class="docutils literal"><span class="pre">solver.u_all</span></tt> and <tt class="docutils literal"><span class="pre">solver.t_all</span></tt>
arrays are available.  These are of interest to see how the adaptive
strategy works between the user-specified time points.</p>
<div class="section" id="the-test-problem">
<h4>The Test Problem<a class="headerlink" href="#the-test-problem" title="Permalink to this headline">¶</a></h4>
<p>We consider the ODE problem for testing adaptive solvers:</p>
<div class="math" id="equation-gaussian:ode:eq">
<span class="eqno">(1)</span>\[     u' = -\frac{t-c}{s^2} u\]</div>
<p>The exact solution is a Gaussian function,</p>
<div class="math">
\[u(t) = \exp{\left(-\frac{1}{2}\left(\frac{t-c}{s}\right)^2\right)}\]</div>
<p>centered around <span class="math">\(t=c\)</span> and width characteristic width (&#8220;standard
deviation&#8221;) <span class="math">\(s\)</span>. The initial condition is taken as the exact <span class="math">\(u\)</span> at <span class="math">\(t=0\)</span>.</p>
<p>Since the Gaussian function is significantly different from zero only in the
interval <span class="math">\([c-3s, c+3s]\)</span>, one may expect that adaptive methods will
efficiently take larger steps when <span class="math">\(u\)</span> is almost constant and increase
the resolution when <span class="math">\(u\)</span> changes substantially in the vicinity
of <span class="math">\(t=c\)</span>. We can test if this is the case with several solvers.</p>
</div>
<div class="section" id="running-simple-methods">
<h4>Running Simple Methods<a class="headerlink" href="#running-simple-methods" title="Permalink to this headline">¶</a></h4>
<p>Let us first use a simple standard method like the 2nd- and 4th-order
Runge-Kutta methods with constant step size. With the former method
(<tt class="docutils literal"><span class="pre">RK2</span></tt>), <span class="math">\(c=3\)</span>, <span class="math">\(s=0.5\)</span>, and <span class="math">\(41\)</span> uniformly distributed points,
the discrepancy between the numerical and exact solution in
Figure <a class="reference internal" href="#gaussian-fig-rk2-41"><em>2nd-order Runge-Kutta method with 41 points for solving :eq:`gaussian:ode:eq`</em></a> is substantial. Increasing the number
of points by a factor of 10 gives a solution much closer to the
exact one, and switching to the 4th-order method (<tt class="docutils literal"><span class="pre">RK4</span></tt>) makes the
curves visually coincide. The problem is therefore quite straightforward
to solve using a sufficient number of points and a higher-order method
(for curiosity we can mention that the <tt class="docutils literal"><span class="pre">ForwardEuler</span></tt> method produces
a maximum value of 0.98 with 20,000 points and 0.998 with 200,000 points).</p>
<div class="figure" id="gaussian-fig-rk2-41">
<img alt="_images/gaussian1_RK2_41.png" src="_images/gaussian1_RK2_41.png" style="width: 500px;" />
<p class="caption"><em>2nd-order Runge-Kutta method with 41 points for solving :eq:`gaussian:ode:eq`</em></p>
</div>
<p>A simple program testing one numerical method goes as follows (<a class="reference external" href="https://github.com/hplgit/odespy/blob/master/doc/src/odespy/src-odespy/gaussian1.py">gaussian1.py</a>).</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">odespy</span><span class="o">,</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">np</span><span class="o">,</span> <span class="nn">matplotlib.pyplot</span> <span class="kn">as</span> <span class="nn">plt</span>

<span class="n">center_point</span> <span class="o">=</span> <span class="mi">3</span>
<span class="n">s</span> <span class="o">=</span> <span class="mf">0.5</span>

<span class="n">problem</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">problems</span><span class="o">.</span><span class="n">Gaussian0</span><span class="p">(</span><span class="n">c</span><span class="o">=</span><span class="n">center_point</span><span class="p">,</span> <span class="n">s</span><span class="o">=</span><span class="n">s</span><span class="p">)</span>

<span class="n">npoints</span> <span class="o">=</span> <span class="mi">41</span>
<span class="n">tp</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">2</span><span class="o">*</span><span class="n">center_point</span><span class="p">,</span> <span class="n">npoints</span><span class="p">)</span>

<span class="n">method</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">RK4</span>
<span class="n">solver</span> <span class="o">=</span> <span class="n">method</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">f</span><span class="p">)</span>
<span class="n">solver</span><span class="o">.</span><span class="n">set_initial_condition</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">U0</span><span class="p">)</span>

<span class="n">u</span><span class="p">,</span> <span class="n">t</span> <span class="o">=</span> <span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">tp</span><span class="p">)</span>

<span class="n">method</span> <span class="o">=</span> <span class="n">solver</span><span class="o">.</span><span class="n">__class__</span><span class="o">.</span><span class="n">__name__</span>
<span class="k">print</span> <span class="s">&#39;</span><span class="si">%.4f</span><span class="s">  </span><span class="si">%s</span><span class="s">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">u</span><span class="o">.</span><span class="n">max</span><span class="p">(),</span> <span class="n">method</span><span class="p">)</span>

<span class="k">if</span> <span class="n">solver</span><span class="o">.</span><span class="n">has_u_t_all</span><span class="p">():</span>
    <span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">solver</span><span class="o">.</span><span class="n">t_all</span><span class="p">,</span> <span class="n">solver</span><span class="o">.</span><span class="n">u_all</span><span class="p">,</span> <span class="s">&#39;bo&#39;</span><span class="p">,</span>
             <span class="n">tp</span><span class="p">,</span> <span class="n">problem</span><span class="o">.</span><span class="n">u_exact</span><span class="p">(</span><span class="n">tp</span><span class="p">))</span>
    <span class="k">print</span> <span class="s">&#39;</span><span class="si">%s</span><span class="s"> used </span><span class="si">%d</span><span class="s"> steps (</span><span class="si">%d</span><span class="s"> specified)&#39;</span> <span class="o">%</span> \
          <span class="p">(</span><span class="n">method</span><span class="p">,</span> <span class="nb">len</span><span class="p">(</span><span class="n">solver</span><span class="o">.</span><span class="n">u_all</span><span class="p">),</span> <span class="nb">len</span><span class="p">(</span><span class="n">tp</span><span class="p">))</span>
<span class="k">else</span><span class="p">:</span>
    <span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">tp</span><span class="p">,</span> <span class="n">solver</span><span class="o">.</span><span class="n">u</span><span class="p">,</span> <span class="n">tp</span><span class="p">,</span> <span class="n">problem</span><span class="o">.</span><span class="n">u_exact</span><span class="p">(</span><span class="n">tp</span><span class="p">))</span>
<span class="n">plt</span><span class="o">.</span><span class="n">legend</span><span class="p">([</span><span class="n">method</span><span class="p">,</span> <span class="s">&#39;exact&#39;</span><span class="p">])</span>
<span class="n">plt</span><span class="o">.</span><span class="n">savefig</span><span class="p">(</span><span class="s">&#39;tmp.png&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="section" id="running-the-runge-kutta-fehlberg-method">
<h4>Running the Runge-Kutta-Fehlberg Method<a class="headerlink" href="#running-the-runge-kutta-fehlberg-method" title="Permalink to this headline">¶</a></h4>
<p>One of the most widely used general-purpose, adaptive methods for
ODE problems is the
Runge-Kutta-Fehlberg method of order (4,5). This method is available
in three alternative implementations in Odespy: a direct Python
version (<tt class="docutils literal"><span class="pre">RKFehlberg</span></tt>), a specialization of a generalized
implementation of explicit adaptive Runge-Kutta methods
(<tt class="docutils literal"><span class="pre">Fehlberg</span></tt>), and as a FORTRAN code (<tt class="docutils literal"><span class="pre">RKF45</span></tt>). We can try one of
these,</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">method</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">Fehlberg</span>
</pre></div>
</div>
<p>Figure <a class="reference internal" href="#gaussian-fig-fehlberg-41"><em>Adaptive Runge-Kutta-Fehlberg method with 41 points for solving :eq:`gaussian:ode:eq`</em></a> shows how <tt class="docutils literal"><span class="pre">Fehlberg</span></tt> with 40 intervals produces a solution of reasonable accuracy. The dots show the actual computational points used by the algorithm (113 adaptively selected points in time).</p>
<div class="figure" id="gaussian-fig-fehlberg-41">
<img alt="_images/gaussian1_Fehlberg_41.png" src="_images/gaussian1_Fehlberg_41.png" style="width: 500px;" />
<p class="caption"><em>Adaptive Runge-Kutta-Fehlberg method with 41 points for solving :eq:`gaussian:ode:eq`</em></p>
</div>
<p>Adaptive algorithms apply an error estimate based on considering a higher-order method as exact, in this case a method of order 5, and a method of lower order (here 4) as the numerically predicted solution. The user can specify an error tolerance. In the program above we just relied to the default tolerance, which can be printed by</p>
<div class="highlight-python"><div class="highlight"><pre><span class="k">print</span> <span class="n">solver</span><span class="o">.</span><span class="n">get</span><span class="p">()</span>
</pre></div>
</div>
<p>yielding a list of all optional parameters:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="p">{</span><span class="s">&#39;f_kwargs&#39;</span><span class="p">:</span> <span class="p">{},</span> <span class="s">&#39;f_args&#39;</span><span class="p">:</span> <span class="p">(),</span>
 <span class="s">&#39;max_step&#39;</span><span class="p">:</span> <span class="mf">1.5000000000000036</span><span class="p">,</span> <span class="s">&#39;verbose&#39;</span><span class="p">:</span> <span class="mi">0</span><span class="p">,</span>
 <span class="s">&#39;min_step&#39;</span><span class="p">:</span> <span class="mf">0.0014999999999999946</span><span class="p">,</span>
 <span class="s">&#39;first_step&#39;</span><span class="p">:</span> <span class="mf">0.14999999999999999</span><span class="p">,</span>
 <span class="s">&#39;rtol&#39;</span><span class="p">:</span> <span class="mf">1e-06</span><span class="p">,</span> <span class="s">&#39;atol&#39;</span><span class="p">:</span> <span class="mf">1e-08</span><span class="p">,</span>
 <span class="s">&#39;name of f&#39;</span><span class="p">:</span> <span class="s">&#39;f&#39;</span><span class="p">,</span> <span class="s">&#39;complex_valued&#39;</span><span class="p">:</span> <span class="bp">False</span><span class="p">,</span>
 <span class="s">&#39;disk_storage&#39;</span><span class="p">:</span> <span class="bp">False</span><span class="p">,</span> <span class="s">&#39;u_exact&#39;</span><span class="p">:</span> <span class="bp">None</span><span class="p">}</span>
</pre></div>
</div>
<p>The tolerances involved are of relative and absolute type, i.e.,</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">estimated_error</span> <span class="o">&lt;=</span> <span class="n">tol</span> <span class="o">=</span> <span class="n">rtol</span><span class="o">*</span><span class="nb">abs</span><span class="p">(</span><span class="n">u</span><span class="p">)</span> <span class="o">+</span> <span class="n">atol</span>
</pre></div>
</div>
<p>is the typical test if the solution is accurate enough.
For very small <tt class="docutils literal"><span class="pre">u</span></tt>, <tt class="docutils literal"><span class="pre">atol</span></tt> comes into play, while for large <tt class="docutils literal"><span class="pre">u</span></tt>, the
relative tolerance <tt class="docutils literal"><span class="pre">rtol</span></tt> dominates.</p>
<p>In this particular example, running <tt class="docutils literal"><span class="pre">RK4</span></tt> with 113 equally spaced points
yields a maximum value of 0.9946, while <tt class="docutils literal"><span class="pre">Fehlberg</span></tt> results in 0.9849.
That is, the much simpler and faster <tt class="docutils literal"><span class="pre">RK4</span></tt> method is also more accurate
than the all-round, &#8220;reliable&#8221; Runge-Kutta-Fehlberg method with an
error tolerance of <span class="math">\(10^-6\)</span>. As we see, the actually error is of the
order <span class="math">\(10^{-2}\)</span> for the latter method.</p>
<p>We can specify stricter tolerances and also control the minimum
allowed step size, <tt class="docutils literal"><span class="pre">min_step</span></tt>, which might be too large to achieve
the desired error level (<a class="reference external" href="https://github.com/hplgit/odespy/blob/master/doc/src/odespy/src-odespy/gaussian2.py">gaussian2.py</a>):</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">rtol</span> <span class="o">=</span> <span class="mf">1E-6</span>
<span class="n">atol</span> <span class="o">=</span> <span class="n">rtol</span>
<span class="n">min_step</span> <span class="o">=</span> <span class="mf">0.000001</span>

<span class="n">solver</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">Fehlberg</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">f</span><span class="p">,</span> <span class="n">atol</span><span class="o">=</span><span class="n">atol</span><span class="p">,</span> <span class="n">rtol</span><span class="o">=</span><span class="n">rtol</span><span class="p">,</span>
                         <span class="n">min_step</span><span class="o">=</span><span class="n">min_step</span><span class="p">)</span>
</pre></div>
</div>
<p>This adjustment does not improve the accuracy. Setting <tt class="docutils literal"><span class="pre">rtol=1E-12</span></tt> does:
the <tt class="docutils literal"><span class="pre">Fehlberg</span></tt> solver then applies 800 points and achieves a maximum
value of 1.00004.</p>
</div>
</div>
<div class="section" id="example-solving-a-complex-ode-problem">
<h3>Example: Solving a Complex ODE Problem<a class="headerlink" href="#example-solving-a-complex-ode-problem" title="Permalink to this headline">¶</a></h3>
<p>Many of the solvers offered by Odespy can deal with complex-valued
ODE problems. Consider</p>
<div class="math">
\[u' = iwu,\quad u(0)=1,\]</div>
<p>where <span class="math">\(i=\sqrt{-1}\)</span> is the imaginary unit.
The right-hand side is implemented as <tt class="docutils literal"><span class="pre">1j*w*u</span></tt> in Python since
Python applies <tt class="docutils literal"><span class="pre">j</span></tt> as the imaginary unit in complex numbers.</p>
<div class="section" id="quick-implementation">
<h4>Quick Implementation<a class="headerlink" href="#quick-implementation" title="Permalink to this headline">¶</a></h4>
<p>For complex-valued ODEs, i.e., complex-valued right-hand side functions
or initial conditions, the argument <tt class="docutils literal"><span class="pre">complex_valued=True</span></tt> must be
supplied to the constructor. A complete program reads</p>
<div class="highlight-python"><div class="highlight"><pre><span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
    <span class="k">return</span> <span class="mi">1</span><span class="n">j</span><span class="o">*</span><span class="n">w</span><span class="o">*</span><span class="n">u</span>

<span class="kn">import</span> <span class="nn">odespy</span><span class="o">,</span> <span class="nn">numpy</span>

<span class="n">w</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="n">numpy</span><span class="o">.</span><span class="n">pi</span>
<span class="n">solver</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">RK4</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="n">complex_valued</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>
<span class="n">solver</span><span class="o">.</span><span class="n">set_initial_condition</span><span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="mi">0</span><span class="n">j</span><span class="p">)</span>
<span class="n">u</span><span class="p">,</span> <span class="n">t</span> <span class="o">=</span> <span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">numpy</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">6</span><span class="p">,</span> <span class="mi">101</span><span class="p">))</span>
</pre></div>
</div>
<p>The function <tt class="docutils literal"><span class="pre">odespy.list_not_suitable_complex_solvers()</span></tt>
returns a list of all the classes in Odespy that are not suitable
for complex-valued ODE problems.</p>
</div>
</div>
<div class="section" id="avoiding-callbacks-to-python">
<h3>Avoiding Callbacks to Python<a class="headerlink" href="#avoiding-callbacks-to-python" title="Permalink to this headline">¶</a></h3>
<p>The ODE solvers that are implemented in FORTRAN calls, by default,
the user&#8217;s Python implementation of <span class="math">\(f(u,t)\)</span>. Making many calls from
FORTRAN to Python may introduce significant overhead and slow down the
solution process. When the algorithm is implemented in FORTRAN we should
also implement the right-hand side in FORTRAN and call this right-hand
side subroutine directly. Odespy offers this possibility.</p>
<p>The idea is that the user writes a FORTRAN subroutine defining <span class="math">\(f(u,t)\)</span>.
Thereafter, <tt class="docutils literal"><span class="pre">f2py</span></tt> is used to make this subroutine callable from Python.
If we specify the Python interface to this subroutine as an <tt class="docutils literal"><span class="pre">f_f77</span></tt>
argument to the solver&#8217;s constructor, the Odespy class will make sure
that no callbacks to the <span class="math">\(f(u,t)\)</span> definition go via Python.</p>
<div class="section" id="the-logistic-ode">
<h4>The Logistic ODE<a class="headerlink" href="#the-logistic-ode" title="Permalink to this headline">¶</a></h4>
<p>Here is a minimalistic example involving the logistic ODE from
the section <a class="reference internal" href="#ode-sec-exgr"><em>First Example: Logistic Growth</em></a>. The FORTRAN implementation of <span class="math">\(f(u,t)\)</span>
is more complicated than the Python counterpart. The subroutine
has the signature</p>
<div class="highlight-fortran"><div class="highlight"><pre>      <span class="k">subroutine </span><span class="nv">f_f77</span><span class="p">(</span><span class="nv">neq</span><span class="p">,</span> <span class="nv">t</span><span class="p">,</span> <span class="nv">u</span><span class="p">,</span> <span class="nv">udot</span><span class="p">)</span>
<span class="nv">Cf2py</span> <span class="k">intent</span><span class="p">(</span><span class="nv">hide</span><span class="p">)</span> <span class="nv">neq</span>
<span class="nv">Cf2py</span> <span class="k">intent</span><span class="p">(</span><span class="nv">out</span><span class="p">)</span> <span class="nv">udot</span>
      <span class="kt">integer </span><span class="nv">neq</span>
      <span class="kt">double precision </span><span class="nv">t</span><span class="p">,</span> <span class="nv">u</span><span class="p">,</span> <span class="nv">udot</span>
      <span class="k">dimension </span><span class="nv">u</span><span class="p">(</span><span class="nv">neq</span><span class="p">),</span> <span class="nv">udot</span><span class="p">(</span><span class="nv">neq</span><span class="p">)</span>
</pre></div>
</div>
<p>This means that there are two additional arguments: <tt class="docutils literal"><span class="pre">neq</span></tt> for the number
of equations in the ODE system, and <tt class="docutils literal"><span class="pre">udot</span></tt> for the array of <span class="math">\(f(u,t)\)</span>
that is output from the subroutine.</p>
<p>We write the FORTRAN implementation of <span class="math">\(f(u,t)\)</span> in a string:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">a</span> <span class="o">=</span> <span class="mi">2</span>
<span class="n">R</span> <span class="o">=</span> <span class="mf">1E+5</span>

<span class="n">f_f77_str</span> <span class="o">=</span> <span class="s">&quot;&quot;&quot;</span>
<span class="s">      subroutine f_f77(neq, t, u, udot)</span>
<span class="s">Cf2py intent(hide) neq</span>
<span class="s">Cf2py intent(out) udot</span>
<span class="s">      integer neq</span>
<span class="s">      double precision t, u, udot</span>
<span class="s">      dimension u(neq), udot(neq)</span>
<span class="s">      udot(1) = </span><span class="si">%.3f</span><span class="s">*u(1)*(1 - u(1)/</span><span class="si">%.1f</span><span class="s">)</span>
<span class="s">      return</span>
<span class="s">      end</span>
<span class="s">&quot;&quot;&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">a</span><span class="p">,</span> <span class="n">R</span><span class="p">)</span>
</pre></div>
</div>
<p>Observe that we can transfer problem parameters to the FORTRAN subroutine
by writing their values directly into the FORTRAN source code.
The other alternative would be to transfer the parameters as global
(COMMON block) variables to the FORTRAN code, which is technically
much more complicated.
Also observe that we need to deal with <tt class="docutils literal"><span class="pre">udot</span></tt> and <tt class="docutils literal"><span class="pre">u</span></tt> as arrays even
for a scalar ODE.</p>
<p>Using <tt class="docutils literal"><span class="pre">f2py</span></tt> to compile the string into a Python module is automated
by the <tt class="docutils literal"><span class="pre">odespy.compile_f77</span></tt> function:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">odespy</span>
<span class="n">f_f77</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">compile_f77</span><span class="p">(</span><span class="n">f_f77_str</span><span class="p">)</span>
</pre></div>
</div>
<p>The returned object <tt class="docutils literal"><span class="pre">f_f77</span></tt> is a callable object that allows the
FORTRAN subroutine to be called as <tt class="docutils literal"><span class="pre">udot</span> <span class="pre">=</span> <span class="pre">f_f77(t,</span> <span class="pre">u)</span></tt> from Python.
(However, the Odespy solvers will not use <tt class="docutils literal"><span class="pre">f_f77</span></tt> directly, but rather
its function pointer to the FORTRAN subroutine, and transfer this pointer
to the FORTRAN solver. The switching between <tt class="docutils literal"><span class="pre">t,</span> <span class="pre">u</span></tt> and <tt class="docutils literal"><span class="pre">u,</span> <span class="pre">t</span></tt> arguments
is taken care of. All necessary steps are automatically done behind
the scene.)</p>
<p>The solver can be declared as</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">solver</span> <span class="o">=</span> <span class="n">odespy</span><span class="o">.</span><span class="n">Lsode</span><span class="p">(</span><span class="n">f</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">f_f77</span><span class="o">=</span><span class="n">f_f77</span><span class="p">)</span>
</pre></div>
</div>
<p>Several solvers accept FORTRAN definitions of the right-hand side:
<tt class="docutils literal"><span class="pre">Lsode</span></tt>, <tt class="docutils literal"><span class="pre">Lsoda</span></tt>, and the other ODEPACK solvers, <tt class="docutils literal"><span class="pre">RKC</span></tt>, <tt class="docutils literal"><span class="pre">RKF45</span></tt>,
<tt class="docutils literal"><span class="pre">Radau5</span></tt>. Look up the documentation of their <tt class="docutils literal"><span class="pre">f_f77</span></tt> parameter to
see exactly what arguments and conventions that the FORTRAN subroutine
demand.</p>
<p>The file <a class="reference external" href="https://github.com/hplgit/odespy/blob/master/doc/src/odespy/src-odespy/logistic10.py">logistic10.py</a> contains a complete program for solving the logistic ODE
with <span class="math">\(f(u,t)\)</span> implemented in Fortran.</p>
</div>
</div>
</div>
</div>


          </div>
        </div>
      </div>
      <div class="sphinxsidebar">
        <div class="sphinxsidebarwrapper">
  <h3><a href="index.html">Table Of Contents</a></h3>
  <ul>
<li><a class="reference internal" href="#"><tt class="docutils literal"><span class="pre">odespy</span></tt></a><ul>
<li><a class="reference internal" href="#basic-usage">Basic Usage</a><ul>
<li><a class="reference internal" href="#overview">Overview</a><ul>
<li><a class="reference internal" href="#step-1">Step 1</a></li>
<li><a class="reference internal" href="#step-2">Step 2</a></li>
<li><a class="reference internal" href="#step-3">Step 3</a></li>
<li><a class="reference internal" href="#step-4">Step 4</a></li>
<li><a class="reference internal" href="#step-5">Step 5</a></li>
<li><a class="reference internal" href="#step-6">Step 6</a></li>
</ul>
</li>
<li><a class="reference internal" href="#first-example-logistic-growth">First Example: Logistic Growth</a></li>
<li><a class="reference internal" href="#parameters-in-the-right-hand-side-function">Parameters in the Right-Hand Side Function</a></li>
<li><a class="reference internal" href="#termination-criterion-for-the-simulation">Termination Criterion for the Simulation</a></li>
<li><a class="reference internal" href="#a-class-based-implementation">A Class-Based Implementation</a></li>
<li><a class="reference internal" href="#using-other-symbols">Using Other Symbols</a></li>
<li><a class="reference internal" href="#example-solving-an-ode-system">Example: Solving an ODE System</a></li>
<li><a class="reference internal" href="#testing-several-methods">Testing Several Methods</a></li>
<li><a class="reference internal" href="#make-a-subclass-of-problem">Make a Subclass of Problem</a></li>
<li><a class="reference internal" href="#using-adaptive-methods">Using Adaptive Methods</a><ul>
<li><a class="reference internal" href="#the-test-problem">The Test Problem</a></li>
<li><a class="reference internal" href="#running-simple-methods">Running Simple Methods</a></li>
<li><a class="reference internal" href="#running-the-runge-kutta-fehlberg-method">Running the Runge-Kutta-Fehlberg Method</a></li>
</ul>
</li>
<li><a class="reference internal" href="#example-solving-a-complex-ode-problem">Example: Solving a Complex ODE Problem</a><ul>
<li><a class="reference internal" href="#quick-implementation">Quick Implementation</a></li>
</ul>
</li>
<li><a class="reference internal" href="#avoiding-callbacks-to-python">Avoiding Callbacks to Python</a><ul>
<li><a class="reference internal" href="#the-logistic-ode">The Logistic ODE</a></li>
</ul>
</li>
</ul>
</li>
</ul>
</li>
</ul>

  <h4>Previous topic</h4>
  <p class="topless"><a href="index.html"
                        title="previous chapter">Odespy API</a></p>
  <h4>Next topic</h4>
  <p class="topless"><a href="rkc.html"
                        title="next chapter"><tt class="docutils literal"><span class="pre">odespy.rkc</span></tt></a></p>
<div id="searchbox" style="display: none">
  <h3>Quick search</h3>
    <form class="search" action="search.html" method="get">
      <input type="text" name="q" />
      <input type="submit" value="Go" />
      <input type="hidden" name="check_keywords" value="yes" />
      <input type="hidden" name="area" value="default" />
    </form>
    <p class="searchtip" style="font-size: 90%">
    Enter search terms or a module, class or function name.
    </p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             >index</a></li>
        <li class="right" >
          <a href="py-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li class="right" >
          <a href="np-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li class="right" >
          <a href="rkc.html" title="odespy.rkc"
             >next</a> |</li>
        <li class="right" >
          <a href="index.html" title="Odespy API"
             >previous</a> |</li>
        <li><a href="index.html">Odespy API 0.2 documentation</a> &raquo;</li> 
      </ul>
    </div>
    <div class="footer">
        &copy; Copyright 2012, Liwei Wang and Hans Petter Langtangen.
      Created using <a href="http://sphinx.pocoo.org/">Sphinx</a> 1.1.2.
    </div>
  </body>
</html>